library(sandwich)
library(mgcv)
library(reshape)
library(plyr)
library(survey)
load("dem.panel.RData")
source("panel-utils.R")


base.vars <- paste(c("camp.length","deminc",
                     "base.poll*deminc","as.factor(year)",
                     "base.und*deminc", "office"
), collapse = "+")
treat.hist <- paste(c("d.gone.neg.l1","d.gone.neg.l2",
                      #"neg.dem.l1", "neg.dem.l2",
                      "d.neg.frac.l3","week*camp.length"#,"as.factor(year)"
),
                    collapse = "+")
covs <- paste(c("dem.polls.l1",
                "neg.rep.l1", "neg.rep.l2", "r.neg.frac.l3",
                "num.rep.l1",
                "num.dem.l1",
                "undother.l1"
), collapse = "+")
cov.names <- c(dem.polls.l1 = "Polls$\\,_{t-1}$",
               neg.rep.l1 = "R Neg$\\,_{t-1}$",
               neg.rep.l2 = "R Neg$\\,_{t-2}$",
               r.neg.frac.l3 = "R Neg Total$\\,_{t-3}$",
               r.neg.frac.l2 = "R Neg Total$\\,_{t-2}$",
               num.rep.l1 = "R Ads$\\,_{t-1}$",
               num.dem.l1 = "D Ads$\\,_{t-1}$",
               undother.l1 = "Undecided$\\,_{t-1}$"
)
n.cov.hist <- paste(c("s(dem.polls.l1, k = 3)",
                      "neg.rep.l1",
                      "neg.rep.l2", "r.neg.frac.l3",
                      "num.rep.l1",
                      "num.dem.l1",
                      "undother.l1"), collapse = "+")


i.cov.hist <- paste(c("neg.rep.l1",
                      #"I(sqrt(num.rep.l1))", #"undother.l1",
                      "neg.rep.l2", #"num.rep.l2",
                      "r.neg.frac.l2",
                      "s(dem.polls.l1)"
), collapse = "+")

i.denom <- as.formula(paste("d.gone.neg ~ ",
                            paste(treat.hist, base.vars,
                                  i.cov.hist, sep = "+")))
n.denom <- as.formula(paste("d.gone.neg ~ ",
                            paste(treat.hist, base.vars,
                                  n.cov.hist, sep = "+")))

numer <- as.formula(paste("d.gone.neg ~ ",
                          paste(treat.hist, base.vars,
                                sep = "+")))

## non.comp are the extremely uncompetitive races that fail common
## support from baseline. Thus, this is an analysis is for the
## population of reasonably competive races.

non.comp <- unique(subset(dem.panel, base.poll > 70 | base.poll < 30
                          | is.na(base.poll) ,
                          "race"))[,1]
dem.panel <- subset(dem.panel, !(race %in% non.comp))

## Positive probability subset remove the weeks where there is either
## no chance of staying positive (too many ads to not have negative
## ones), or no chance of going negative (race is too safe, or there
## are no ads)

dem.panel$pos.prob.subset <- (dem.panel$week !=0 &
                                dem.panel$week > dem.panel$first.week +2 &
                                dem.panel$num.dem != 0 )




dem.panel$dem.pscore <- NA
dem.panel$dem.marpscore <- NA
dem.panel$sw <- NA

dem.wmod.ni <- iptw(denominator = n.denom, numerator = numer,
                    data = dem.panel, id = "race", time = "week",
                    family = "binomial", subset = deminc == 0,
                    ones  = !dem.panel$pos.prob.subset,
                    bal.covs = strsplit(covs, "\\+")[[1]])

dem.wmod.inc <- iptw(denominator = i.denom, numerator = numer,
                     data = dem.panel, id = "race", time = "week",
                     family = "binomial", subset = deminc == 1,
                     ones  = !dem.panel$pos.prob.subset,
                     bal.covs = strsplit(covs, "\\+")[[1]])
summary(dem.wmod.ni)
summary(dem.wmod.inc)
dem.panel$sw[dem.panel$deminc == 1] <- dem.wmod.inc$sw
dem.panel$sw[dem.panel$deminc == 0] <- dem.wmod.ni$sw
dem.panel$dem.pscore[dem.panel$deminc == 1] <- dem.wmod.inc$d.pscore
dem.panel$dem.pscore[dem.panel$deminc == 0] <- dem.wmod.ni$d.pscore


dem.final <- subset(dem.panel, week == -1)

causal.des <- svydesign(ids = ~ race, weights = ~ sw, data =
                          dem.final[!is.na(dem.final$sw),])


share.form <- paste("turnout ~ ",
                    "deminc * d.neg.rec + deminc * I(d.neg.dur - d.neg.rec)+as.factor(year)+",
                    base.vars)

share.form <- as.formula(share.form)


causal.share <-svyglm(share.form, causal.des)
summary(causal.share)



inc.form <- paste("turnout ~ ",
                  "I(deminc == 0) * d.neg.rec   + I(deminc == 0) * I(d.neg.dur - d.neg.rec) +",
                  gsub("deminc", "I(deminc == 0)", base.vars))

inc.form <- as.formula(inc.form)
causal.share.inc <- svyglm(inc.form, causal.des)
summary(causal.share.inc)



sims <- 500
confound <- seq(-4, 4, by = 0.2)
sens.sha <- matrix(NA, nrow = length(confound), ncol =  sims)
sens.inc <- matrix(NA, nrow = length(confound), ncol = sims)
sens.early.sha <- sens.sha
sens.early.inc <- sens.inc


strategic.bias <- function(alpha, a) return(alpha)
confirmation.bias <- function(alpha, a) {alpha * (2 * a - 1)}


dem.panel$y.a <- NULL
for (i in 1:length(confound)) {
  cat(i, "\n")
  this.alpha <- confound[i]
  dem.panel <- ddply(dem.panel, "race", transform,
                     tmp = sens.adj(y=turnout, a=d.gone.neg,
                                    a.fits=dem.pscore, alpha=this.alpha, confound.func = confirmation.bias))
  newNames <- names(dem.panel)
  names(newNames) <- names(dem.panel)
  newNames["tmp"] <- paste("y.a.", i, sep = "")
  names(dem.panel) <- newNames
}


dem.final <- dem.panel[dem.panel$week == -1, ]
dem.final <- dem.final[!is.na(dem.final$sw),]




race.list <- unique(dem.final$race[!is.na(dem.final$sw)])
race.lengths <- c()
race.rows <- list()
for (i in 1: length(race.list)) {
  race.rows[[i]] <- rownames(dem.panel[which(dem.panel$race
                                             ==
                                               race.list[i]),])
  race.lengths <- c(race.lengths, nrow(dem.panel[which(dem.panel$race ==
                                                         race.list[i]),]))

}
names(race.rows) <- race.list
names(race.lengths) <- race.list

for (i in 1:sims) {
  cat(i, " ")
  if (i %% 10 == 0) cat("\n")
  races <- sample(race.list, replace = TRUE,
                  size = length(race.list))
  races <- as.character(races)
  star <- unlist(race.rows[races])
  bootid <- rep(1:length(races), times = race.lengths[races])
  dem.star <- dem.panel[star,]
  dem.star$bootid <- bootid

  dem.wmod.ni <- iptw(denominator = n.denom, numerator = numer,
                      data = dem.star, id = "bootid", time = "week",
                      family = "binomial", subset = deminc == 0,
                      ones  = !dem.star$pos.prob.subset,
                      bal.covs = strsplit(covs, "\\+")[[1]])

  dem.wmod.inc <- iptw(denominator = i.denom, numerator = numer,
                       data = dem.star, id = "bootid", time = "week",
                       family = "binomial", subset = deminc == 1,
                       ones  = !dem.star$pos.prob.subset,
                       bal.covs = strsplit(covs, "\\+")[[1]])
  dem.star$sw[dem.star$deminc == 1] <- dem.wmod.inc$sw
  dem.star$sw[dem.star$deminc == 0] <- dem.wmod.ni$sw
  dem.star$dem.pscore[dem.star$deminc == 1] <- dem.wmod.inc$d.pscore
  dem.star$dem.pscore[dem.star$deminc == 0] <- dem.wmod.ni$d.pscore

  dem.star <- dem.star[dem.star$week == -1,]
  dem.star <- dem.star[!is.na(dem.star$sw),]

  causal.des.star <- svydesign(ids = ~ race, weights = ~ sw,
                               data = dem.star)


  for (ii in 1:length(confound)) {

    share.form <- paste("y.a.", ii, " ~ ",
                        "deminc * d.neg.rec + deminc * I(d.neg.dur - d.neg.rec) + ",
                        base.vars, sep = "")

    share.form <- as.formula(share.form)


    inc.form <- paste("y.a.", ii, " ~ ",
                      "I(deminc == 0) * d.neg.rec + I(deminc == 0) * I(d.neg.dur - d.neg.rec) +  ",
                      base.vars, sep = "")

    inc.form <- as.formula(inc.form)


    causal.sha.star <- svyglm(share.form, design = causal.des.star)
    causal.inc.star <- svyglm(inc.form, design = causal.des.star)



    sens.sha[ii, i] <- causal.sha.star$coef["d.neg.rec"]
    sens.inc[ii, i] <- causal.inc.star$coef["d.neg.rec"]
    sens.early.sha[ii,i] <- causal.sha.star$coef["I(d.neg.dur - d.neg.rec)"]
    sens.early.inc[ii,i] <-
      causal.inc.star$coef["I(d.neg.dur - d.neg.rec)"]
  }


}

causal.des <- svydesign(ids = ~ race, weights = ~ sw, data =
                          dem.final[!is.na(dem.final$sw),])

rsqs <- rep(NA, length(confound))
for (ii in 1:length(confound)) {
  rsqs.form <- paste("I(y.a.", ii,")",
                     " ~  deminc * I(d.neg.dur - d.neg.rec) + ",
                      base.vars, sep = "")

  rsqs.form <- as.formula(rsqs.form)


  tmp <- svyglm(rsqs.form, design = causal.des)
  rsqs[ii] <- confound[ii]^2*var(dem.final$d.neg.rec)/(var(tmp$residuals))
}


##save(list = c("rsqs", "sens.sha", "sens.inc", "confound"),
##     file = "turnout-dyn.RData")

##load("turnout-dyn.RData")

regTable <- each(mean, sd, lower.95.bound, upper.95.bound)
noninc <- adply(sens.sha, 1, regTable)[,-1]
inc <- adply(sens.inc, 1, regTable)[,-1]
colnames(noninc) <- c("Mean", "SE", "95% Lower Bound",
                      "95% Upper Bound")
colnames(inc) <- c("Mean", "SE", "95% Lower Bound",
                   "95% Upper Bound")

noninc$pos <- noninc[,3] > 0 & noninc[,4] > 0
noninc$non <- noninc[,3] < 0 & noninc[,4] > 0
noninc$neg <- noninc[,3] < 0 & noninc[,4] < 0
noninc.pos <- c(which(noninc$pos), max(which(noninc$pos)) + 1)
noninc.neg <- c(min(which(noninc$neg))-1, which(noninc$neg))
inc$pos <- inc[,3] > 0 & inc[,4] > 0
inc$non <- inc[,3] < 0 & inc[,4] > 0
inc$neg <- inc[,3] < 0 & inc[,4] < 0
inc.pos <- c(which(inc$pos), max(which(inc$pos)) + 1)
inc.neg <- c(min(which(inc$neg))-1, which(inc$neg))

rsqs <- sign(confound) * rsqs

### Sensitivity Figure - raw scale ###
##tikz(file = "../figs/turnout-dyn-raw.tex", width = 5.5, height = 3)
par(mfrow = c(1,2), cex.main = .72, cex.axis = 0.72, cex.lab = 0.72,
    mar = c(4.2,4,1,1) + .1)
plot(x = NULL, xlim = c(-4, 4), ylim = c(-4, 5), axes = FALSE,
     xlab = "Amount of confounding (α)", main = "Democrat Nonincumbents",
     ylab = "Estimated effect")

axis(side = 1, ps = 4)
axis(side = 2, las = 2)
abline(v = 0, col = rgb(0.1,0.1,0.1, alpha = 0.25))
polygon(x = c(confound[noninc.pos],rev(confound[noninc.pos])), y = c(noninc[noninc.pos, 3], rev(noninc[noninc.pos, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                                    alpha = 0.25), border =
          FALSE)
polygon(x = c(confound[noninc$non],rev(confound[noninc$non])), y = c(noninc[noninc$non, 3], rev(noninc[noninc$non, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                                    alpha = 0.25), border =
          FALSE, density = 25)
polygon(x = c(confound[noninc.neg],rev(confound[noninc.neg])), y = c(noninc[noninc.neg, 3], rev(noninc[noninc.neg, 4])) , col = rgb(0.5, 0.5, 0.5,
                                                                                                                                    alpha = 0.25), border =
          FALSE)
lines(x = confound, y = noninc[,1], col = rgb(0.1,0.1,0.1),
      lwd = 2, lty = 2)
abline(h = 0, col = rgb(0.1, 0.1, 0.1, alpha = .25), lwd = 1)

text(x = 0, y = -4, "Negative campaigns\n lower turnout", pos = 2,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))

text(x = 0, y = -4, "Negative campaigns\n higher turnout", pos = 4,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))
arrows(x0 = c(-0.2, 0.2), x1 = c(-2, 2), y0 = c(-3, -3), y1 =
         c(-3,-3), col = rgb(0.5, 0.5, 0.5, alpha = .7), length =
         0.05)


plot(x = NULL, xlim = c(-4, 4), ylim = c(-4, 5), axes = FALSE,
     xlab = "Amount of confounding (α)", main = "Democrat Incumbents",
     ylab = "Estimated effect", cex = .72)
axis(side = 1)
axis(side = 2, las = 2)
abline(v = 0, col = rgb(0.1,0.1,0.1, alpha = 0.25))
polygon(x = c(confound[inc.pos],rev(confound[inc.pos])), y = c(inc[inc.pos, 3], rev(inc[inc.pos, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                  alpha = 0.25), border =
          FALSE)
polygon(x = c(confound[36:41],rev(confound[36:41])), y = c(inc[36:41, 3], rev(inc[36:41, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                      alpha = 0.25), border =
          FALSE, density = 25)
## polygon(x = c(confound[29:41],rev(confound[29:41])), y = c(inc[29:41, 3], rev(inc[29:41, 4])) , col = rgb(0.05, 0.05, 0.05,
##                                                                                                   alpha = 0.25), border =
##         FALSE, density = 25)

polygon(x = c(confound[inc.neg],rev(confound[inc.neg])), y = c(inc[inc.neg, 3], rev(inc[inc.neg, 4])) , col = rgb(0.5, 0.5, 0.5,
                                                                                                                  alpha = 0.25), border =
          FALSE)
lines(x = confound, y = inc[,1], col = rgb(0.1,0.1,0.1),
      lwd = 2, lty = 2)
abline(h = 0, col = rgb(0.1, 0.1, 0.1, alpha = .25), lwd = 1)

text(x = 0, y = -4, "Negative campaigns\n lower turnout", pos = 2,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))

text(x = 0, y = -4, "Negative campaigns\n higher turnout", pos = 4,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))
arrows(x0 = c(-0.2, 0.2), x1 = c(-2, 2), y0 = c(-3, -3), y1 =
         c(-3,-3), col = rgb(0.5, 0.5, 0.5, alpha = .7), length =
         0.05)

##dev.off()

### Sensitivity Figure - R^2 scale ###
##tikz(file = "../figs/turnout-dyn.tex", width = 5.5, height = 3)
par(mfrow = c(1,2), cex.main = .8, cex.axis = 1, cex.lab = 0.8,
    mar = c(4.2,4,1,1) + .1)
plot(x = NULL, xlim = c(-.75, .75), ylim = c(-4, 5), axes = FALSE,
     xlab = "Variance explained by confounding \n($R_α^2$)", main = "Democrat Nonincumbents",
     ylab = "Estimated effect")

axis(side = 1, ps = 4)
axis(side = 2, las = 2)
abline(v = 0, col = rgb(0.1,0.1,0.1, alpha = 0.25))
polygon(x = c(rsqs[noninc.pos],rev(rsqs[noninc.pos])), y = c(noninc[noninc.pos, 3], rev(noninc[noninc.pos, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                                    alpha = 0.25), border =
          FALSE)
polygon(x = c(rsqs[noninc$non],rev(rsqs[noninc$non])), y = c(noninc[noninc$non, 3], rev(noninc[noninc$non, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                                    alpha = 0.25), border =
          FALSE, density = 25)
polygon(x = c(rsqs[noninc.neg],rev(rsqs[noninc.neg])), y = c(noninc[noninc.neg, 3], rev(noninc[noninc.neg, 4])) , col = rgb(0.5, 0.5, 0.5,
                                                                                                                                    alpha = 0.25), border =
          FALSE)
lines(x = rsqs, y = noninc[,1], col = rgb(0.1,0.1,0.1),
      lwd = 2, lty = 2)
abline(h = 0, col = rgb(0.1, 0.1, 0.1, alpha = .25), lwd = 1)

text(x = 0, y = -4, "Negative campaigns\n lower turnout", pos = 2,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))

text(x = 0, y = -4, "Negative campaigns\n higher turnout", pos = 4,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))
arrows(x0 = c(-0.1, 0.1), x1 = c(-.4, 0.4), y0 = c(-2.8,-2.8), y1 =
         c(-2.8,-2.8), col = rgb(0.5, 0.5, 0.5, alpha = .7), length =
         0.05)


plot(x = NULL, xlim = c(-.75, 0.75), ylim = c(-4, 5), axes = FALSE,
     xlab = "Variance explained by confounding \n($R_α^2$)", main = "Democrat Incumbents",
     ylab = "Estimated effect", cex = .72)
axis(side = 1)
axis(side = 2, las = 2)
abline(v = 0, col = rgb(0.1,0.1,0.1, alpha = 0.25))
polygon(x = c(rsqs[inc.pos],rev(rsqs[inc.pos])), y = c(inc[inc.pos, 3], rev(inc[inc.pos, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                                  alpha = 0.25), border =
          FALSE)
polygon(x = c(rsqs[36:41],rev(rsqs[36:41])), y = c(inc[36:41, 3], rev(inc[36:41, 4])) , col = rgb(0.05, 0.05, 0.05,
                                                                                                          alpha = 0.25), border =
          FALSE, density = 25)
## polygon(x = c(rsqs[29:41],rev(rsqs[29:41])), y = c(inc[29:41, 3], rev(inc[29:41, 4])) , col = rgb(0.05, 0.05, 0.05,
##                                                                                                   alpha = 0.25), border =
##         FALSE, density = 25)

polygon(x = c(rsqs[inc.neg],rev(rsqs[inc.neg])), y = c(inc[inc.neg, 3], rev(inc[inc.neg, 4])) , col = rgb(0.5, 0.5, 0.5,
                                                                                                                  alpha = 0.25), border =
          FALSE)
lines(x = rsqs, y = inc[,1], col = rgb(0.1,0.1,0.1),
      lwd = 2, lty = 2)
abline(h = 0, col = rgb(0.1, 0.1, 0.1, alpha = .25), lwd = 1)

text(x = 0, y = -4, "Negative campaigns\n lower turnout", pos = 2,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))

text(x = 0, y = -4, "Negative campaigns\n higher turnout", pos = 4,
     cex = 0.7, col = rgb(0.5, 0.5, 0.5))
arrows(x0 = c(-0.1, 0.1), x1 = c(-0.4, 0.4), y0 = c(-2.8, -2.8), y1 =
         c(-2.8,-2.8), col = rgb(0.5, 0.5, 0.5, alpha = .7), length =
         0.05)

##dev.off()
